R Stan examples: IQ recovery after TBI
Georges Monette
York UniversityMarch 29, 2017 at 13:47
Stan for Asymptotic Recovery Curves
We can use non-linear mixed effects models for normally distributed responses to estimate asymptotic recovery curves for VIQ and PIQ separately. It’s feasible to extend these models to handle VIQ and PIQ as a multivariate response, but doing so is not straighforward.
Here, we use Stan to build a model for a multivariate response which will allow us to compare the parameters for the two models in a way that exploits the covariance between the responses.
We start with univariate models which we then combine into a multivariate model.
Loading packages:
library(spida2)
library(magrittr, pos = 1000)
library(lattice)
library(latticeExtra)
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())Explore data
data(iq)
?iq
head(iq) DAYSPC DCOMA SEX AgeIn ID PIQ VIQ
1 30 4 Male 20.67077 3358 87 89
2 16 17 Male 55.28816 3535 95 77
3 40 1 Male 55.91513 3547 95 116
4 13 10 Male 61.66461 3592 59 73
5 19 6 Male 30.12731 3728 67 73
6 13 3 Male 57.06229 3790 76 69names(iq) <- tolower(names(iq))
dd <-iq
(p <- xyplot(piq ~ dayspc | sex, dd, groups = id, type = 'b'))
update(p, xlim = c(0,4000))
if(interactive) {
rows <- numeric(0)
# can repeat:
trellis.focus()
rows <- c(rows, panel.identify(labels=dd$id))
# end
trellis.unfocus()
rows
save(rows, file = 'rows.rda')
}
load('rows.rda', verbose = T)Loading objects:
rowsdd[rows,] %>% sortdf(~dcoma+dayspc) dayspc dcoma sex agein id piq viq
47 3333 9 Male 43.93977 2600 86 80
327 3337 9 Male 43.93977 2600 101 84
303 1491 21 Male 22.00684 651 71 94
325 3412 21 Male 22.00684 651 68 92
310 1926 130 Male 28.27379 1939 95 108
321 3111 130 Male 28.27379 1939 88 111
326 3864 130 Male 28.27379 1939 88 105# id = 2600 retested 4 days apart
dd[dd$id %in% dd[rows,'id'],] %>% sortdf(~dcoma+dayspc) dayspc dcoma sex agein id piq viq
47 3333 9 Male 43.93977 2600 86 80
327 3337 9 Male 43.93977 2600 101 84
303 1491 21 Male 22.00684 651 71 94
325 3412 21 Male 22.00684 651 68 92
234 295 130 Male 28.27379 1939 67 117
271 562 130 Male 28.27379 1939 85 111
310 1926 130 Male 28.27379 1939 95 108
321 3111 130 Male 28.27379 1939 88 111
326 3864 130 Male 28.27379 1939 88 105if(interactive){
library(p3d)
Init3d()
dd$dcoma.cat <- cut(dd$dcoma, c(-1,2,5,10,20,50,Inf))
Plot3d( viq ~ piq + log(dayspc) |
dcoma.cat, dd, groups = id,
col = heat.colors(6))
Plot3d( viq ~ log(dcoma+2) + log(dayspc) |
dcoma.cat, dd, groups = id,
col = heat.colors(6))
Plot3d( piq ~ log(dcoma+2) + log(dayspc) |
dcoma.cat, dd, groups = id,
col = heat.colors(12)[1:6])
fg()
Id3d()
}Stan program
stanfile <- 'asymp_model_1.stan'
pr <- function(x) print(readLines(x), quote = FALSE) # write a function for repetitive tasks
pr(stanfile) [1]
[2] // Asymptotic recovery model
[3] // - random intercept
[4]
[5] data {
[6] int N;
[7] int J;
[8] vector[N] y;
[9] vector[N] time;
[10] vector[N] coma;
[11] int id[N];
[12] int Np; // predictor values for prediction
[13] vector[Np] time_p;
[14] vector[Np] coma_p;
[15] }
[16] transformed data{
[17] real ln2;
[18] ln2 = log(2); // factor for half-recovery time
[19] }
[20] parameters {
[21] real hrt; // half-recovery time
[22] real asymp; // mean asymptotic recovery if dcoma = 0
[23] real bcoma; // duration of coma parameter
[24] real init_def; // initial deficit
[25] vector[J] u; // random intercept residual
[26] real <lower=0> sigma; // sd within
[27] real <lower=0> sigma_u; // sd between
[28] }
[29] transformed parameters {
[30] real estd_reliability;
[31] estd_reliability = sigma_u^2 / (sigma_u^2 + sigma^2) ;
[32] }
[33] model {
[34] u ~ normal(0,sigma_u);
[35] y ~ normal(asymp + u[id] + bcoma * coma +
[36] init_def * exp(-time/(hrt*ln2)), sigma);
[37] }
[38] generated quantities {
[39] vector[Np] y_fit;
[40] y_fit = asymp + bcoma * coma_p +
[41] init_def * exp(-time_p/(hrt*ln2));
[42] } system.time(
asymp_model_dso <- stan_model(stanfile)
) user system elapsed
0.094 0.035 0.842 Prepare data list
names(dd)[1] "dayspc" "dcoma" "sex" "agein" "id" "piq" "viq" pred <- expand.grid(time=seq(0, 3*365, 30),
coma = seq(0,30,5))
dat <- list(
N = nrow(dd),
id = nid <- as.numeric(as.factor(dd$id)),
J = max(nid),
y = dd$piq,
time = dd$dayspc,
coma = sqrt(dd$dcoma),
Np = nrow(pred),
time_p = pred$time,
coma_p = sqrt(pred$coma)
)
set.seed(789723)
mod <- sampling(asymp_model_dso, dat, chains = 4, iter = 2000)
# Note: chains and iter are defaults
# traceplot(mod)
names(mod) [1] "hrt" "asymp" "bcoma"
[4] "init_def" "u[1]" "u[2]"
[7] "u[3]" "u[4]" "u[5]"
[10] "u[6]" "u[7]" "u[8]"
[13] "u[9]" "u[10]" "u[11]"
[16] "u[12]" "u[13]" "u[14]"
[19] "u[15]" "u[16]" "u[17]"
[22] "u[18]" "u[19]" "u[20]"
[25] "u[21]" "u[22]" "u[23]"
[28] "u[24]" "u[25]" "u[26]"
[31] "u[27]" "u[28]" "u[29]"
[34] "u[30]" "u[31]" "u[32]"
[37] "u[33]" "u[34]" "u[35]"
[40] "u[36]" "u[37]" "u[38]"
[43] "u[39]" "u[40]" "u[41]"
[46] "u[42]" "u[43]" "u[44]"
[49] "u[45]" "u[46]" "u[47]"
[52] "u[48]" "u[49]" "u[50]"
[55] "u[51]" "u[52]" "u[53]"
[58] "u[54]" "u[55]" "u[56]"
[61] "u[57]" "u[58]" "u[59]"
[64] "u[60]" "u[61]" "u[62]"
[67] "u[63]" "u[64]" "u[65]"
[70] "u[66]" "u[67]" "u[68]"
[73] "u[69]" "u[70]" "u[71]"
[76] "u[72]" "u[73]" "u[74]"
[79] "u[75]" "u[76]" "u[77]"
[82] "u[78]" "u[79]" "u[80]"
[85] "u[81]" "u[82]" "u[83]"
[88] "u[84]" "u[85]" "u[86]"
[91] "u[87]" "u[88]" "u[89]"
[94] "u[90]" "u[91]" "u[92]"
[97] "u[93]" "u[94]" "u[95]"
[100] "u[96]" "u[97]" "u[98]"
[103] "u[99]" "u[100]" "u[101]"
[106] "u[102]" "u[103]" "u[104]"
[109] "u[105]" "u[106]" "u[107]"
[112] "u[108]" "u[109]" "u[110]"
[115] "u[111]" "u[112]" "u[113]"
[118] "u[114]" "u[115]" "u[116]"
[121] "u[117]" "u[118]" "u[119]"
[124] "u[120]" "u[121]" "u[122]"
[127] "u[123]" "u[124]" "u[125]"
[130] "u[126]" "u[127]" "u[128]"
[133] "u[129]" "u[130]" "u[131]"
[136] "u[132]" "u[133]" "u[134]"
[139] "u[135]" "u[136]" "u[137]"
[142] "u[138]" "u[139]" "u[140]"
[145] "u[141]" "u[142]" "u[143]"
[148] "u[144]" "u[145]" "u[146]"
[151] "u[147]" "u[148]" "u[149]"
[154] "u[150]" "u[151]" "u[152]"
[157] "u[153]" "u[154]" "u[155]"
[160] "u[156]" "u[157]" "u[158]"
[163] "u[159]" "u[160]" "u[161]"
[166] "u[162]" "u[163]" "u[164]"
[169] "u[165]" "u[166]" "u[167]"
[172] "u[168]" "u[169]" "u[170]"
[175] "u[171]" "u[172]" "u[173]"
[178] "u[174]" "u[175]" "u[176]"
[181] "u[177]" "u[178]" "u[179]"
[184] "u[180]" "u[181]" "u[182]"
[187] "u[183]" "u[184]" "u[185]"
[190] "u[186]" "u[187]" "u[188]"
[193] "u[189]" "u[190]" "u[191]"
[196] "u[192]" "u[193]" "u[194]"
[199] "u[195]" "u[196]" "u[197]"
[202] "u[198]" "u[199]" "u[200]"
[205] "sigma" "sigma_u" "estd_reliability"
[208] "y_fit[1]" "y_fit[2]" "y_fit[3]"
[211] "y_fit[4]" "y_fit[5]" "y_fit[6]"
[214] "y_fit[7]" "y_fit[8]" "y_fit[9]"
[217] "y_fit[10]" "y_fit[11]" "y_fit[12]"
[220] "y_fit[13]" "y_fit[14]" "y_fit[15]"
[223] "y_fit[16]" "y_fit[17]" "y_fit[18]"
[226] "y_fit[19]" "y_fit[20]" "y_fit[21]"
[229] "y_fit[22]" "y_fit[23]" "y_fit[24]"
[232] "y_fit[25]" "y_fit[26]" "y_fit[27]"
[235] "y_fit[28]" "y_fit[29]" "y_fit[30]"
[238] "y_fit[31]" "y_fit[32]" "y_fit[33]"
[241] "y_fit[34]" "y_fit[35]" "y_fit[36]"
[244] "y_fit[37]" "y_fit[38]" "y_fit[39]"
[247] "y_fit[40]" "y_fit[41]" "y_fit[42]"
[250] "y_fit[43]" "y_fit[44]" "y_fit[45]"
[253] "y_fit[46]" "y_fit[47]" "y_fit[48]"
[256] "y_fit[49]" "y_fit[50]" "y_fit[51]"
[259] "y_fit[52]" "y_fit[53]" "y_fit[54]"
[262] "y_fit[55]" "y_fit[56]" "y_fit[57]"
[265] "y_fit[58]" "y_fit[59]" "y_fit[60]"
[268] "y_fit[61]" "y_fit[62]" "y_fit[63]"
[271] "y_fit[64]" "y_fit[65]" "y_fit[66]"
[274] "y_fit[67]" "y_fit[68]" "y_fit[69]"
[277] "y_fit[70]" "y_fit[71]" "y_fit[72]"
[280] "y_fit[73]" "y_fit[74]" "y_fit[75]"
[283] "y_fit[76]" "y_fit[77]" "y_fit[78]"
[286] "y_fit[79]" "y_fit[80]" "y_fit[81]"
[289] "y_fit[82]" "y_fit[83]" "y_fit[84]"
[292] "y_fit[85]" "y_fit[86]" "y_fit[87]"
[295] "y_fit[88]" "y_fit[89]" "y_fit[90]"
[298] "y_fit[91]" "y_fit[92]" "y_fit[93]"
[301] "y_fit[94]" "y_fit[95]" "y_fit[96]"
[304] "y_fit[97]" "y_fit[98]" "y_fit[99]"
[307] "y_fit[100]" "y_fit[101]" "y_fit[102]"
[310] "y_fit[103]" "y_fit[104]" "y_fit[105]"
[313] "y_fit[106]" "y_fit[107]" "y_fit[108]"
[316] "y_fit[109]" "y_fit[110]" "y_fit[111]"
[319] "y_fit[112]" "y_fit[113]" "y_fit[114]"
[322] "y_fit[115]" "y_fit[116]" "y_fit[117]"
[325] "y_fit[118]" "y_fit[119]" "y_fit[120]"
[328] "y_fit[121]" "y_fit[122]" "y_fit[123]"
[331] "y_fit[124]" "y_fit[125]" "y_fit[126]"
[334] "y_fit[127]" "y_fit[128]" "y_fit[129]"
[337] "y_fit[130]" "y_fit[131]" "y_fit[132]"
[340] "y_fit[133]" "y_fit[134]" "y_fit[135]"
[343] "y_fit[136]" "y_fit[137]" "y_fit[138]"
[346] "y_fit[139]" "y_fit[140]" "y_fit[141]"
[349] "y_fit[142]" "y_fit[143]" "y_fit[144]"
[352] "y_fit[145]" "y_fit[146]" "y_fit[147]"
[355] "y_fit[148]" "y_fit[149]" "y_fit[150]"
[358] "y_fit[151]" "y_fit[152]" "y_fit[153]"
[361] "y_fit[154]" "y_fit[155]" "y_fit[156]"
[364] "y_fit[157]" "y_fit[158]" "y_fit[159]"
[367] "y_fit[160]" "y_fit[161]" "y_fit[162]"
[370] "y_fit[163]" "y_fit[164]" "y_fit[165]"
[373] "y_fit[166]" "y_fit[167]" "y_fit[168]"
[376] "y_fit[169]" "y_fit[170]" "y_fit[171]"
[379] "y_fit[172]" "y_fit[173]" "y_fit[174]"
[382] "y_fit[175]" "y_fit[176]" "y_fit[177]"
[385] "y_fit[178]" "y_fit[179]" "y_fit[180]"
[388] "y_fit[181]" "y_fit[182]" "y_fit[183]"
[391] "y_fit[184]" "y_fit[185]" "y_fit[186]"
[394] "y_fit[187]" "y_fit[188]" "y_fit[189]"
[397] "y_fit[190]" "y_fit[191]" "y_fit[192]"
[400] "y_fit[193]" "y_fit[194]" "y_fit[195]"
[403] "y_fit[196]" "y_fit[197]" "y_fit[198]"
[406] "y_fit[199]" "y_fit[200]" "y_fit[201]"
[409] "y_fit[202]" "y_fit[203]" "y_fit[204]"
[412] "y_fit[205]" "y_fit[206]" "y_fit[207]"
[415] "y_fit[208]" "y_fit[209]" "y_fit[210]"
[418] "y_fit[211]" "y_fit[212]" "y_fit[213]"
[421] "y_fit[214]" "y_fit[215]" "y_fit[216]"
[424] "y_fit[217]" "y_fit[218]" "y_fit[219]"
[427] "y_fit[220]" "y_fit[221]" "y_fit[222]"
[430] "y_fit[223]" "y_fit[224]" "y_fit[225]"
[433] "y_fit[226]" "y_fit[227]" "y_fit[228]"
[436] "y_fit[229]" "y_fit[230]" "y_fit[231]"
[439] "y_fit[232]" "y_fit[233]" "y_fit[234]"
[442] "y_fit[235]" "y_fit[236]" "y_fit[237]"
[445] "y_fit[238]" "y_fit[239]" "y_fit[240]"
[448] "y_fit[241]" "y_fit[242]" "y_fit[243]"
[451] "y_fit[244]" "y_fit[245]" "y_fit[246]"
[454] "y_fit[247]" "y_fit[248]" "y_fit[249]"
[457] "y_fit[250]" "y_fit[251]" "y_fit[252]"
[460] "y_fit[253]" "y_fit[254]" "y_fit[255]"
[463] "y_fit[256]" "y_fit[257]" "y_fit[258]"
[466] "y_fit[259]" "lp__" pars <- grepv('^u|y_fit', names(mod), invert = T)
pars[1] "hrt" "asymp" "bcoma"
[4] "init_def" "sigma" "sigma_u"
[7] "estd_reliability" "lp__" traceplot(mod, pars = pars)
pairs(mod, pars = pars)
if(interactive) {
library(shinystan)
mod_sso <- launch_shinystan(mod)
}
zd <- as.data.frame(mod, pars = pars)
head(zd) hrt asymp bcoma init_def sigma sigma_u
1 352.0651 102.14958 -2.409008 -17.19386 7.255366 11.40271
2 280.1005 100.73336 -2.460842 -17.74048 7.549275 14.95228
3 202.0059 96.23615 -1.543932 -18.50024 6.701158 12.54657
4 254.4917 100.87878 -1.832296 -20.22619 6.789999 12.59399
5 284.2113 100.09366 -2.042092 -17.07062 7.634592 13.40865
6 247.9092 96.78811 -1.456427 -17.42886 6.868731 12.95669
estd_reliability lp__
1 0.7118160 -1426.493
2 0.7968666 -1443.540
3 0.7780491 -1416.420
4 0.7747862 -1414.997
5 0.7551779 -1417.888
6 0.7806170 -1401.731xyplot(hrt~asymp, zd)
print(mod, pars = pars)Inference for Stan model: asymp_model_1.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50%
hrt 277.27 1.07 64.93 173.72 231.85 267.65
asymp 100.96 0.05 1.99 97.13 99.63 100.92
bcoma -1.92 0.01 0.40 -2.69 -2.20 -1.93
init_def -19.50 0.04 1.87 -23.25 -20.73 -19.50
sigma 6.92 0.01 0.44 6.13 6.61 6.90
sigma_u 13.19 0.01 0.84 11.59 12.62 13.16
estd_reliability 0.78 0.00 0.03 0.71 0.76 0.78
lp__ -1415.25 0.52 14.77 -1444.76 -1425.33 -1414.83
75% 97.5% n_eff Rhat
hrt 315.79 430.61 3692 1.00
asymp 102.29 104.94 1567 1.00
bcoma -1.64 -1.13 1342 1.00
init_def -18.27 -15.75 2800 1.00
sigma 7.21 7.84 1486 1.00
sigma_u 13.74 14.95 4000 1.00
estd_reliability 0.81 0.84 1834 1.00
lp__ -1404.74 -1388.23 810 1.01
Samples were drawn using NUTS(diag_e) at Wed Mar 29 13:47:50 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).get_posterior_mean(mod, pars = pars)[,5] hrt asymp bcoma init_def
277.267594 100.955620 -1.920296 -19.503431
sigma sigma_u estd_reliability lp__
6.921886 13.193070 0.782237 -1415.249592 Note that you can do extensive plotting of the posterior sample with ‘shinystan’ including some not highly satisfying 3d plots.
Random initial deficit
It wasn’t feasible to fit a model that included random initial deficits using ‘lme’ so this will be an interesting experiment.
Note that the data will be identical.
stanfile <- 'asymp_model_2.stan'
pr(stanfile) [1]
[2] // Asymptotic recovery model
[3] // - random intercept
[4]
[5] data {
[6] int N;
[7] int J;
[8] vector[N] y;
[9] vector[N] time;
[10] vector[N] coma;
[11] int id[N];
[12] int Np; // predictor values for prediction
[13] vector[Np] time_p;
[14] vector[Np] coma_p;
[15] }
[16] transformed data{
[17] real ln2;
[18] ln2 = log(2); // factor for half-recovery time
[19] }
[20] parameters {
[21] real hrt; // half-recovery time
[22] real asymp; // mean asymptotic recovery if dcoma = 0
[23] real bcoma; // duration of coma parameter
[24] real init_def_mean; // initial deficit
[25] vector[J] u; // random intercept residual
[26] vector[J] u_init_def; // random initial deficit
[27] real <lower=0> sigma; // sd within
[28] real <lower=0> sigma_u; // sd between
[29] real <lower=0> sigma_idef;
[30] }
[31] transformed parameters {
[32] real estd_reliability;
[33] estd_reliability = sigma_u^2 / (sigma_u^2 + sigma^2) ;
[34] }
[35] model {
[36] // note few changes from non-random init_def
[37] u ~ normal(0,sigma_u);
[38] // add an empirical prior for u_init_def:
[39] u_init_def ~ normal(init_def_mean, sigma_idef);
[40] // add u_init_def[id] and element-wise multiplication
[41] // of vectors (.*):
[42] y ~ normal(asymp + u[id] + bcoma * coma +
[43] (init_def_mean + u_init_def[id])
[44] .* exp(-time/(hrt*ln2)), sigma);
[45] }
[46] generated quantities {
[47] vector[Np] y_fit;
[48] y_fit = asymp + bcoma * coma_p +
[49] init_def_mean * exp(-time_p/(hrt*ln2));
[50] } system.time(
asymp_model_2_dso <- stan_model(stanfile)
)In file included from /usr/local/lib/R/site-library/BH/include/boost/config.hpp:39:0,
from /usr/local/lib/R/site-library/BH/include/boost/math/tools/config.hpp:13,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/var.hpp:7,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core.hpp:12,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/mat.hpp:4,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math.hpp:4,
from /usr/local/lib/R/site-library/StanHeaders/include/src/stan/model/model_header.hpp:4,
from file385317a558f2.cpp:8:
/usr/local/lib/R/site-library/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined [enabled by default]
# define BOOST_NO_CXX11_RVALUE_REFERENCES
^
<command-line>:0:0: note: this is the location of the previous definition user system elapsed
33.396 2.446 36.522 Fitting the model. The data is the same.
mod2 <- sampling(asymp_model_2_dso, dat, chains = 4, iter = 2000)Warning: There were 745 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupWarning: There were 3 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-lowWarning: Examine the pairs() plot to diagnose sampling problemspars <- grepv('^u|y_fit', names(mod2), invert = T)
pars[1] "hrt" "asymp" "bcoma"
[4] "init_def_mean" "sigma" "sigma_u"
[7] "sigma_idef" "estd_reliability" "lp__" traceplot(mod2, pars = pars)
pairs(mod2, pars = pars)
Poor convergence: possible remedies
The variability in the sd of initial deficit is poorly estimated.
We can try to improve this in a least two ways:
- Reparametrize the model: Taking the initial point at time 0 puts it outside the range of the data. Using a later time as time 0 could improve the estimate. This can be done just by changing the data and subtracting a constant from time. We then need to remember the changed interpretation of parameters although most parameters are not affected, e.g. ‘asymp’ and ‘hrt’
- Use an informative prior for ‘sigma_idef’. So far all models have used flat priors for hyper-parameters. Gelman would not hesitate to use much more informative priors.
We’ll try step 1 because it’s easy. We’ll subtract 180 from time so that time 0 will be relocated to approximately 6 months.
Reparametrize
dat_6 <- list(
N = nrow(dd),
id = nid <- as.numeric(as.factor(dd$id)),
J = max(nid),
y = dd$piq,
time = dd$dayspc - 180,
coma = sqrt(dd$dcoma),
Np = nrow(pred),
time_p = pred$time - 180,
coma_p = sqrt(pred$coma)
)
mod2_6 <- sampling(asymp_model_2_dso, dat_6, chains = 4, iter = 2000)Warning: There were 815 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupWarning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-lowWarning: Examine the pairs() plot to diagnose sampling problemspars <- grepv('^u|y_fit', names(mod2_6), invert = T)
pars[1] "hrt" "asymp" "bcoma"
[4] "init_def_mean" "sigma" "sigma_u"
[7] "sigma_idef" "estd_reliability" "lp__" traceplot(mod2_6, pars = pars)
pairs(mod2_6, pars = pars)
Informative prior
We will try a Gamma prior on ‘sigma_idef’
stanfile <- 'asymp_model_3.stan'
pr(stanfile) [1]
[2] // Asymptotic recovery model
[3] // - random intercept
[4]
[5] data {
[6] int N;
[7] int J;
[8] vector[N] y;
[9] vector[N] time;
[10] vector[N] coma;
[11] int id[N];
[12] int Np; // predictor values for prediction
[13] vector[Np] time_p;
[14] vector[Np] coma_p;
[15] // gamma prior parameters
[16] real gamma_df;
[17] real gamma_scale;
[18] }
[19] transformed data{
[20] real ln2;
[21] ln2 = log(2); // factor for half-recovery time
[22] }
[23] parameters {
[24] real hrt; // half-recovery time
[25] real asymp; // mean asymptotic recovery if dcoma = 0
[26] real bcoma; // duration of coma parameter
[27] real init_def_mean; // initial deficit
[28] vector[J] u; // random intercept residual
[29] vector[J] u_init_def; // random initial deficit
[30] real <lower=0> sigma; // sd within
[31] real <lower=0> sigma_u; // sd between
[32] real <lower=0> sigma_idef;
[33] }
[34] transformed parameters {
[35] real estd_reliability;
[36] estd_reliability = sigma_u^2 / (sigma_u^2 + sigma^2) ;
[37] }
[38] model {
[39] // Gamma prior
[40] sigma_idef ~ gamma(gamma_df, 1/gamma_scale);
[41] // note few changes from non-random init_def
[42] u ~ normal(0,sigma_u);
[43] // add an empirical prior for u_init_def:
[44] u_init_def ~ normal(init_def_mean, sigma_idef);
[45] // add u_init_def[id] and element-wise multiplication
[46] // of vectors (.*):
[47] y ~ normal(asymp + u[id] + bcoma * coma +
[48] (init_def_mean + u_init_def[id])
[49] .* exp(-time/(hrt*ln2)), sigma);
[50] }
[51] generated quantities {
[52] vector[Np] y_fit;
[53] y_fit = asymp + bcoma * coma_p +
[54] init_def_mean * exp(-time_p/(hrt*ln2));
[55] } system.time(
asymp_model_3_dso <- stan_model(stanfile)
)In file included from /usr/local/lib/R/site-library/BH/include/boost/config.hpp:39:0,
from /usr/local/lib/R/site-library/BH/include/boost/math/tools/config.hpp:13,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/var.hpp:7,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core.hpp:12,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/mat.hpp:4,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math.hpp:4,
from /usr/local/lib/R/site-library/StanHeaders/include/src/stan/model/model_header.hpp:4,
from file385352cf9003.cpp:8:
/usr/local/lib/R/site-library/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined [enabled by default]
# define BOOST_NO_CXX11_RVALUE_REFERENCES
^
<command-line>:0:0: note: this is the location of the previous definition user system elapsed
33.660 0.914 35.301 Fitting the model. The data is the same plus the two parameters for the gamma prior.
mod3_6 <- sampling(asymp_model_3_dso,
c(dat_6, gamma_df = 4, gamma_scale = .5),
chains = 4, iter = 2000,
control = list(adapt_delta = 0.85))Warning: There were 87 divergent transitions after warmup. Increasing adapt_delta above 0.85 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupWarning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-lowWarning: Examine the pairs() plot to diagnose sampling problemspars <- grepv('^u|y_fit', names(mod3_6), invert = T)
pars[1] "hrt" "asymp" "bcoma"
[4] "init_def_mean" "sigma" "sigma_u"
[7] "sigma_idef" "estd_reliability" "lp__" traceplot(mod3_6, pars = pars)
pairs(mod3_6, pars = pars)Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
print(mod3_6, pars = pars)Inference for Stan model: asymp_model_3.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5%
hrt 1.56740e+02 1.032100e+02 1.521100e+02 6.750000e+00
asymp 5.09500e+01 3.539000e+01 5.008000e+01 1.700000e-01
bcoma -9.20000e-01 6.900000e-01 1.010000e+00 -2.540000e+00
init_def_mean -2.62000e+00 9.700000e-01 1.460000e+00 -5.300000e+00
sigma 3.75000e+00 2.220000e+00 3.160000e+00 1.400000e-01
sigma_u 8.87000e+00 3.210000e+00 4.570000e+00 2.410000e+00
sigma_idef 1.60000e+00 8.700000e-01 1.280000e+00 2.500000e-01
estd_reliability 8.80000e-01 7.000000e-02 1.000000e-01 7.200000e-01
lp__ -1.38287e+32 1.554271e+32 4.849026e+32 -2.718517e+33
25% 50% 75% 97.5% n_eff Rhat
hrt 1.448000e+01 7.911000e+01 290.33 424.99 2 3.33
asymp 1.260000e+00 4.822000e+01 100.99 104.18 2 38.71
bcoma -1.890000e+00 -2.700000e-01 -0.01 0.14 2 3.68
init_def_mean -3.910000e+00 -1.640000e+00 -1.43 -1.00 2 2.81
sigma 8.700000e-01 3.330000e+00 6.85 7.65 2 10.26
sigma_u 5.720000e+00 8.720000e+00 13.11 14.52 2 8.29
sigma_idef 2.500000e-01 1.240000e+00 3.58 3.58 2 3.73
estd_reliability 7.900000e-01 9.300000e-01 0.98 1.00 2 4.42
lp__ -2.679679e+31 -3.924177e+11 -1556.21 -1362.73 10 1.38
Samples were drawn using NUTS(diag_e) at Wed Mar 29 13:54:15 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).get_posterior_mean(mod3_6, pars = pars)[,5] %>% cbind .
hrt 1.567368e+02
asymp 5.094795e+01
bcoma -9.159970e-01
init_def_mean -2.620795e+00
sigma 3.747419e+00
sigma_u 8.870717e+00
sigma_idef 1.603489e+00
estd_reliability 8.841998e-01
lp__ -1.382870e+32get_posterior_mean(mod2_6, pars = pars)[,5] %>% cbind .
hrt 1.778507e+04
asymp 9.013245e+01
bcoma -1.981526e+00
init_def_mean -1.050906e+01
sigma 6.558797e+00
sigma_u 8.059216e+00
sigma_idef 2.496417e+00
estd_reliability 4.656921e-01
lp__ -3.148898e+31get_posterior_mean(mod2, pars = pars)[,5] %>% cbind .
hrt 276.8901869
asymp 100.6649413
bcoma -1.8387664
init_def_mean -9.7724863
sigma 6.8567078
sigma_u 13.1091332
sigma_idef 3.7960164
estd_reliability 0.7827252
lp__ -1733.0739153Baysian uniform prior
We’ll just use our intuition about IQs. The standard deviation couldn’t ‘possibly’ be less than 3 and couldn’t ‘possibly’ be greater than 30!
stanfile <- 'asymp_model_3b.stan'
pr(stanfile) [1]
[2] // Asymptotic recovery model
[3] // - random intercept
[4] // - 'Bayesian' prior on sigma_idef
[5]
[6] data {
[7] int N;
[8] int J;
[9] vector[N] y;
[10] vector[N] time;
[11] vector[N] coma;
[12] int id[N];
[13] int Np; // predictor values for prediction
[14] vector[Np] time_p;
[15] vector[Np] coma_p;
[16] }
[17] transformed data{
[18] real ln2;
[19] ln2 = log(2); // factor for half-recovery time
[20] }
[21] parameters {
[22] real hrt; // half-recovery time
[23] real asymp; // mean asymptotic recovery if dcoma = 0
[24] real bcoma; // duration of coma parameter
[25] real init_def_mean; // initial deficit
[26] vector[J] u; // random intercept residual
[27] vector[J] u_init_def; // random initial deficit
[28] real <lower=0> sigma; // sd within
[29] real <lower=0> sigma_u; // sd between
[30] real <lower=3,upper=30> sigma_idef; // uniform prior from 3 to 30
[31] }
[32] transformed parameters {
[33] real estd_reliability;
[34] estd_reliability = sigma_u^2 / (sigma_u^2 + sigma^2) ;
[35] }
[36] model {
[37] // Gamma prior
[38] // sigma_idef ~ gamma(gamma_df, 1/gamma_scale); // commenting out the gamma prior
[39] // note few changes from non-random init_def
[40] u ~ normal(0,sigma_u);
[41] // add an empirical prior for u_init_def:
[42] u_init_def ~ normal(init_def_mean, sigma_idef);
[43] // add u_init_def[id] and element-wise multiplication
[44] // of vectors (.*):
[45] y ~ normal(asymp + u[id] + bcoma * coma +
[46] (init_def_mean + u_init_def[id])
[47] .* exp(-time/(hrt*ln2)), sigma);
[48] }
[49] generated quantities {
[50] vector[Np] y_fit;
[51] y_fit = asymp + bcoma * coma_p +
[52] init_def_mean * exp(-time_p/(hrt*ln2));
[53] } system.time(
asymp_model_3b_dso <- stan_model(stanfile)
)In file included from /usr/local/lib/R/site-library/BH/include/boost/config.hpp:39:0,
from /usr/local/lib/R/site-library/BH/include/boost/math/tools/config.hpp:13,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/var.hpp:7,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core.hpp:12,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/mat.hpp:4,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math.hpp:4,
from /usr/local/lib/R/site-library/StanHeaders/include/src/stan/model/model_header.hpp:4,
from file385363a7af4e.cpp:8:
/usr/local/lib/R/site-library/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined [enabled by default]
# define BOOST_NO_CXX11_RVALUE_REFERENCES
^
<command-line>:0:0: note: this is the location of the previous definition user system elapsed
33.790 0.897 35.409 Fitting the model. The data is the same plus the two parameters for the gamma prior.
mod3b_6 <- sampling(asymp_model_3b_dso,
dat_6,
chains = 4, iter = 2000,
control = list(adapt_delta = 0.85))Warning: There were 691 divergent transitions after warmup. Increasing adapt_delta above 0.85 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupWarning: There were 3 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-lowWarning: Examine the pairs() plot to diagnose sampling problemspars <- grepv('^u|y_fit', names(mod3b_6), invert = T)
pars[1] "hrt" "asymp" "bcoma"
[4] "init_def_mean" "sigma" "sigma_u"
[7] "sigma_idef" "estd_reliability" "lp__" traceplot(mod3b_6, pars = pars)
pairs(mod3b_6, pars = pars)Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
print(mod3b_6, pars = pars)Inference for Stan model: asymp_model_3b.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5%
hrt 3.356300e+02 2.641500e+02 6.938500e+02 2.590000e+00
asymp 5.138000e+01 3.605000e+01 5.102000e+01 1.500000e-01
bcoma -1.280000e+00 6.700000e-01 9.800000e-01 -2.520000e+00
init_def_mean -1.960000e+00 2.590000e+00 3.720000e+00 -7.610000e+00
sigma 1.889000e+01 1.668000e+01 2.768000e+01 4.900000e-01
sigma_u 7.170000e+00 3.550000e+00 5.070000e+00 1.400000e-01
sigma_idef 8.900000e+00 3.460000e+00 4.990000e+00 3.050000e+00
estd_reliability 6.000000e-01 2.600000e-01 3.700000e-01 0.000000e+00
lp__ -6.936327e+80 7.594519e+80 1.612219e+81 -6.538357e+81
25% 50% 75% 97.5% n_eff Rhat
hrt 4.153000e+01 137.78 522.59 1038.61 7 1.22
asymp 5.200000e-01 48.31 102.34 106.14 2 37.69
bcoma -1.920000e+00 -1.75 -0.14 0.43 2 3.83
init_def_mean -5.520000e+00 -0.68 1.68 1.68 2 5.84
sigma 5.000000e-01 7.14 8.66 79.36 3 2.74
sigma_u 3.270000e+00 4.37 12.10 14.20 2 7.41
sigma_idef 3.900000e+00 9.40 13.29 15.81 2 7.76
estd_reliability 1.400000e-01 0.73 0.91 0.99 2 9.36
lp__ -1.156439e+80 -2078.99 -1795.19 -1710.55 5 2.32
Samples were drawn using NUTS(diag_e) at Wed Mar 29 13:56:50 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).get_posterior_mean(mod3b_6, pars = pars)[,5] %>% cbind .
hrt 3.356334e+02
asymp 5.137587e+01
bcoma -1.284750e+00
init_def_mean -1.959387e+00
sigma 1.889447e+01
sigma_u 7.166087e+00
sigma_idef 8.904122e+00
estd_reliability 6.036895e-01
lp__ -6.936327e+80get_posterior_mean(mod2_6, pars = pars)[,5] %>% cbind .
hrt 1.778507e+04
asymp 9.013245e+01
bcoma -1.981526e+00
init_def_mean -1.050906e+01
sigma 6.558797e+00
sigma_u 8.059216e+00
sigma_idef 2.496417e+00
estd_reliability 4.656921e-01
lp__ -3.148898e+31get_posterior_mean(mod2, pars = pars)[,5] %>% cbind .
hrt 276.8901869
asymp 100.6649413
bcoma -1.8387664
init_def_mean -9.7724863
sigma 6.8567078
sigma_u 13.1091332
sigma_idef 3.7960164
estd_reliability 0.7827252
lp__ -1733.0739153Alas, this didn’t work either.
Sensitivity analysis
Since it looks like we can’t estimate ‘sigma_idef’ from the data with our model but we are concerned that it might have an uncertain impact on our other results, we can resort to a sensitivity analysis. Try a variety of plausible values and see how it affects other results.
We also put a plausible lower bound of 10 on sigma_u.
stanfile <- 'asymp_model_3c.stan'
pr(stanfile) [1]
[2] // Asymptotic recovery model
[3] // - random intercept
[4] // - sensitivity analysis for sigma_idef
[5]
[6] data {
[7] int N;
[8] int J;
[9] vector[N] y;
[10] vector[N] time;
[11] vector[N] coma;
[12] int <lower=1,upper=J> id[N];
[13] int Np; // predictor values for prediction
[14] vector[Np] time_p;
[15] vector[Np] coma_p;
[16] real sigma_idef;
[17] }
[18] transformed data{
[19] real ln2;
[20] ln2 = log(2); // factor for half-recovery time
[21] }
[22] parameters {
[23] real hrt; // half-recovery time
[24] real asymp; // mean asymptotic recovery if dcoma = 0
[25] real bcoma; // duration of coma parameter
[26] real init_def_mean; // initial deficit
[27] vector[J] u; // random intercept residual
[28] vector[J] u_init_def; // random initial deficit
[29] real <lower=0.0> sigma; // sd within
[30] real <lower=10.0> sigma_u; // sd between // lower limit on sigma_u
[31] // real <lower=3,upper=30> sigma_idef; // will input as data for sensitivity analysis
[32] }
[33] transformed parameters {
[34] real estd_reliability;
[35] estd_reliability = sigma_u^2 / (sigma_u^2 + sigma^2) ;
[36] }
[37] model {
[38] // Gamma prior
[39] // sigma_idef ~ gamma(gamma_df, 1/gamma_scale); // commenting out the gamma prior
[40] // note few changes from non-random init_def
[41] u ~ normal(0,sigma_u);
[42] // add an empirical prior for u_init_def:
[43] u_init_def ~ normal(init_def_mean, sigma_idef);
[44] // add u_init_def[id] and element-wise multiplication
[45] // of vectors (.*):
[46] y ~ normal(asymp + u[id] + bcoma * coma +
[47] (init_def_mean + u_init_def[id])
[48] .* exp(-time/(hrt*ln2)), sigma);
[49] }
[50] generated quantities {
[51] vector[Np] y_fit;
[52] y_fit = asymp + bcoma * coma_p +
[53] init_def_mean * exp(-time_p/(hrt*ln2));
[54] } system.time(
asymp_model_3c_dso <- stan_model(stanfile)
) user system elapsed
0.098 0.030 0.837 Fitting the model. The data is the same plus the two parameters for the gamma prior.
mod3c_6 <- sampling(asymp_model_3c_dso,
c(dat_6, sigma_idef = .1),
chains = 4, iter = 2000)Warning: There were 72 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupWarning: There were 1 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-lowWarning: Examine the pairs() plot to diagnose sampling problemspars <- grepv('^u|y_fit', names(mod3c_6), invert = T)
pars[1] "hrt" "asymp" "bcoma"
[4] "init_def_mean" "sigma" "sigma_u"
[7] "estd_reliability" "lp__" traceplot(mod3c_6, pars = pars)
pairs(mod3c_6, pars = pars)
print(mod3c_6, pars = pars)Inference for Stan model: asymp_model_3c.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5%
hrt 2.186900e+02 8.997000e+01 1.358900e+02 2.740000e+00
asymp 7.633000e+01 3.039000e+01 4.301000e+01 1.900000e+00
bcoma -1.680000e+00 3.400000e-01 5.900000e-01 -2.730000e+00
init_def_mean -2.560000e+00 1.610000e+00 2.340000e+00 -5.140000e+00
sigma 5.320000e+00 1.980000e+00 2.830000e+00 4.700000e-01
sigma_u 1.301000e+01 2.100000e-01 7.600000e-01 1.173000e+01
estd_reliability 8.400000e-01 7.000000e-02 1.000000e-01 7.100000e-01
lp__ -4.106383e+76 4.787314e+76 8.285904e+76 -2.918382e+77
25% 50% 75% 97.5% n_eff Rhat
hrt 1.020900e+02 259.77 313.00 417.63 2 2.68
asymp 7.105000e+01 100.31 101.92 104.76 2 27.26
bcoma -2.130000e+00 -1.75 -0.89 -0.89 3 1.67
init_def_mean -4.220000e+00 -3.56 -0.73 1.33 2 3.89
sigma 4.380000e+00 6.72 7.12 7.81 2 7.69
sigma_u 1.254000e+01 12.81 13.51 14.73 13 1.08
estd_reliability 7.700000e-01 0.80 0.91 1.00 2 3.46
lp__ -1.159598e+76 -1524.22 -1508.49 -1484.01 3 3.38
Samples were drawn using NUTS(diag_e) at Wed Mar 29 13:59:04 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).get_posterior_mean(mod3c_6, pars = pars)[,5] %>% cbind .
hrt 2.186874e+02
asymp 7.633317e+01
bcoma -1.680734e+00
init_def_mean -2.564241e+00
sigma 5.323000e+00
sigma_u 1.301404e+01
estd_reliability 8.352306e-01
lp__ -4.106383e+76get_posterior_mean(mod3c_6, pars = pars)[,5] %>% cbind .
hrt 2.186874e+02
asymp 7.633317e+01
bcoma -1.680734e+00
init_def_mean -2.564241e+00
sigma 5.323000e+00
sigma_u 1.301404e+01
estd_reliability 8.352306e-01
lp__ -4.106383e+76Multivariate model for VIQ and PIQ
Instead of simple standard deviations for each VIQ and PIQ, we need covariance matrices to describe the within- and the between-subject covariance between VIQ and PIQ.
The first model uses a uniform prior on covariance matrices using constraints generated by Stan on covariance matrices.
It is more efficient to use a combination of an LKJ prior on correlation together with priors of variances. The second model below use this form of parametrization.
stanfile <- 'asymp_model_4.stan'
pr(stanfile) [1] //
[2] // Multivariate model for VIQ and PIQ
[3] //
[4]
[5] data {
[6] int N;
[7] int J;
[8] matrix[N,2] iq;
[9] vector[N] time;
[10] vector[N] coma;
[11] int id[N];
[12] }
[13] transformed data{
[14] real ln2;
[15] vector[2] zero;
[16] ln2 = log(2);
[17] for ( i in 1:2) zero[i] = 0.0;
[18] }
[19] parameters {
[20] vector <lower=1,upper=10000>[2] hrt;
[21] vector <lower=0,upper=200>[2] asymp;
[22] vector <lower=-100,upper=100>[2] init_def;
[23] vector [2] bcoma;
[24] vector[2] u[J];
[25] cov_matrix[2] Sigma;
[26] cov_matrix[2] Sigma_u;
[27] }
[28] transformed parameters {
[29] real hrt_diff;
[30] real bcoma_diff;
[31] hrt_diff = hrt[2] - hrt[1];
[32] bcoma_diff = bcoma[2] - bcoma[1];
[33] }
[34] model {
[35] vector[2] eta;
[36] // for the multinormal distribution we need to loop over observations
[37] for(j in 1:J) u[j] ~ multi_normal(zero, Sigma_u);
[38] for(n in 1:N) {
[39] eta[1] = asymp[1] + u[id[n],1] + bcoma[1] * coma[n] +
[40] init_def[1] * exp(-time[n]/(hrt[1]*ln2));
[41] eta[2] = asymp[2] + u[id[n],2] + bcoma[2] * coma[n] +
[42] init_def[2] * exp(-time[n]/(hrt[2]*ln2));
[43] iq[n,] ~ multi_normal(eta, Sigma);
[44] }
[45] }
[46] system.time(
asymp_model_4_dso <- stan_model(stanfile)
)In file included from /usr/local/lib/R/site-library/BH/include/boost/config.hpp:39:0,
from /usr/local/lib/R/site-library/BH/include/boost/math/tools/config.hpp:13,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/var.hpp:7,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core.hpp:12,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/mat.hpp:4,
from /usr/local/lib/R/site-library/StanHeaders/include/stan/math.hpp:4,
from /usr/local/lib/R/site-library/StanHeaders/include/src/stan/model/model_header.hpp:4,
from file38533f792847.cpp:8:
/usr/local/lib/R/site-library/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined [enabled by default]
# define BOOST_NO_CXX11_RVALUE_REFERENCES
^
<command-line>:0:0: note: this is the location of the previous definition user system elapsed
36.302 1.233 38.270 Data for multivariate model
dat <- list(
N = nrow(dd),
id = nid <- as.numeric(as.factor(dd$id)),
J = max(nid),
iq = cbind(dd$viq,dd$piq),
time = dd$dayspc,
coma = sqrt(dd$dcoma)
)
mult_mod <- sampling( asymp_model_4_dso, dat)
pars <- grepv('^u' ,names(mult_mod), invert = T)
traceplot(mult_mod, pars = pars)
pairs(mult_mod, pars = pars)the following parameters were dropped because they are duplicative
Sigma[1,2] Sigma_u[1,2]
print(mult_mod, pars = pars)Inference for Stan model: asymp_model_4.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75%
hrt[1] 66.00 0.45 19.35 36.99 52.62 63.39 76.24
hrt[2] 246.89 0.94 51.88 161.55 209.98 242.37 278.12
asymp[1] 99.69 0.06 1.60 96.70 98.58 99.65 100.75
asymp[2] 100.35 0.07 1.94 96.58 99.06 100.33 101.64
init_def[1] -23.23 0.11 4.74 -34.41 -25.74 -22.62 -19.94
init_def[2] -19.43 0.05 2.02 -23.14 -20.78 -19.43 -18.12
bcoma[1] -0.68 0.02 0.39 -1.45 -0.94 -0.68 -0.43
bcoma[2] -1.91 0.02 0.42 -2.74 -2.19 -1.91 -1.63
Sigma[1,1] 33.10 0.11 4.35 25.73 29.97 32.70 35.88
Sigma[2,1] 20.33 0.12 4.29 12.73 17.39 20.02 22.80
Sigma[1,2] 20.33 0.12 4.29 12.73 17.39 20.02 22.80
Sigma[2,2] 49.60 0.18 6.66 38.24 45.06 48.87 53.62
Sigma_u[1,1] 162.71 0.39 19.54 128.64 149.15 161.34 175.20
Sigma_u[2,1] 119.30 0.38 17.84 87.32 107.00 118.07 130.66
Sigma_u[1,2] 119.30 0.38 17.84 87.32 107.00 118.07 130.66
Sigma_u[2,2] 175.87 0.46 22.68 135.33 160.06 174.46 190.14
hrt_diff 180.89 0.80 50.81 98.17 145.43 175.23 210.50
bcoma_diff -1.22 0.01 0.33 -1.85 -1.45 -1.23 -1.00
lp__ -2603.56 0.71 21.04 -2645.33 -2618.05 -2602.95 -2589.04
97.5% n_eff Rhat
hrt[1] 111.44 1864 1.00
hrt[2] 362.33 3063 1.00
asymp[1] 102.93 620 1.01
asymp[2] 104.26 804 1.00
init_def[1] -15.79 1844 1.00
init_def[2] -15.30 1658 1.00
bcoma[1] 0.11 657 1.01
bcoma[2] -1.08 774 1.00
Sigma[1,1] 42.72 1495 1.00
Sigma[2,1] 29.84 1256 1.00
Sigma[1,2] 29.84 1256 1.00
Sigma[2,2] 64.49 1394 1.00
Sigma_u[1,1] 203.72 2518 1.00
Sigma_u[2,1] 156.86 2226 1.00
Sigma_u[1,2] 156.86 2226 1.00
Sigma_u[2,2] 224.79 2440 1.00
hrt_diff 295.17 4000 1.00
bcoma_diff -0.58 3180 1.00
lp__ -2564.57 877 1.00
Samples were drawn using NUTS(diag_e) at Wed Mar 29 14:05:08 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).